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^ Abstract. We apply a variational technique to solve the time-dependent Gross- 

' Pitaevskii equation for Bose-Einstein condensates in which an additional dipole-dipole 

I— i, interaction between the atoms is present with the goal of modelling the dynamics of 

i-^ ' such condensates. We show that universal stability thresholds for the collapse of the 

condensates correspond to bifurcation points where always two stationary solutions 
pH , of the Gross-Pitaevskii equation disappear in a tangent bifurcation, one dynamically 

d ' stable and the other unstable. We point out that the thresholds also correspond to 

^ ^. "exceptional points" , i.e. branching singularities of the Hamiltonian. We analyse the 

1 1 1 dynamics of excited condensate wave functions via Poincare surfaces of section for the 

condensate parameters and find both regular and chaotic motion, corresponding to 
(quasi-) periodically oscillating and irregularly fluctuating condensates, respectively. 

1/^ ' Stable islands are found to persist up to energies well above the saddle point of the 

. mean- field energy, alongside with collapsing modes. The results are applicable when 

I the shape of the condensate is axisymmetric. 
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■ At sufficiently low temperatures a condensate of weakly interacting bosons can 

^ ■ be represented by a single wave function whose dynamics obeys the dynamics of the 
■ - - ' Gross-Pitaevskii equation [1,2]. The equation is nonlinear in the wave function, and 
therefore the solutions of the equation exhibit features not familiar from solutions of 
ordinary Schrodinger equations of quantum mechanics. As an example of the effects 
of the nonlinearity, Huepe et al. [3,4] demonstrated that for Bose-Einstein condensates 
with attractive contact interaction, described by a negative s-wave scattering length a, 
bifurcations of the stationary solutions of the Gross-Pitaevskii equation appear. These 
authors also determined both the stable (elliptic) and the unstable (hyperbolic) branches 
of the solutions. In physical terms, the bifurcation points correspond to critical particle 
numbers, above which, for given strength of the attractive interaction, collapse of the 
condensate sets in. For Bose-Einstein condensates of ^Li [5, 6] and ^^Rb atoms [7, 8] 
these collapses were experimentally observed. 

In those condensates the short-range contact interaction is the only interaction to be 
considered. By contrast, in Bose-Einstein condensates of dipolar gases [9-13] also a long- 
range dipole-dipole interaction is present. This offers the unique opportunity to study 
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degenerate quantum gases with adjustable long-range and short-range interactions. The 
achievement of Bose-Einstein condensation in a gas of chromium atoms [14], with a 
large dipole moment, has in fact opened the way to promising experiments on dipolar 
gases [15], which could show a wealth of novel phenomena [16-19]. In particular, 
the experimental observation of the collapse of dipolar quantum gases has also been 
reported [20] which occurs when the contact interaction is reduced, for a given particle 
number, below some critical value using a Feshbach resonance. 

In this experimental situation it is most timely and appropriate to extend the 
investigations of the effects of the nonlinearity of the Gross-Pitaevskii equation to 
dipolar quantum gases, and this is the goal of the present paper. To model these 
effects we will pursue, for the sake of simplicity, a variational ansatz. We do this in 
the spirit of Refs. [3, 4, 21] where also variational techniques were applied to model 
the dynamics of dilute ultracold atom clouds in the Bose-Einstein condensed phase by 
solving the Gross-Pitaevskii equation without dipole-dipole interaction. In fact, quite 
recently Parker et al. [22] have pointed out that in dipolar Bose-Einstein condensates the 
Gaussian variational method gives excellent agreement with the full numerical solutions 
of the Gross-Pitaevskii equation in wide ranges of the physical parameters. For oblate 
dipolar Bose-Einstein condensates the Gaussian approximation appears to agree only 
for weak dipolar interactions. The approximation used may not be valid in the limit 
of interest and clearly further exact studies based on exact solutions of the dipolar 
Gross-Pitaevskii equation need to be carried out to verify the results. Full numerical 
quantum calculations, for condensates with only the contact interaction [3,4,21] present, 
or with an additional attractive gravity-like 1/r interaction [23-26], have confirmed that 
properties of the solutions of the Gross-Pitaevskii equation found in the variational 
calculations are recovered in the quantum calculations. We should therefore expect 
that a simple variational approach will also capture essential features of the dynamics 
of condensate wave functions of dipolar gases. 

We treat the problem in the "atomic" units provided by the magnetic dipole-dipole 
interaction, i.e., we measure lengths in units of the "dipole length" 



energies in units of -Ed = ^^/(2?Tiaj), frequencies in units of cUd = E^^/h and time in 
units of h/E^, respectively. In ([T]), a is the fine-structure constant, Oq the Bohr radius, 
m/me the ratio of the atom and electron mass, and /i the magnetic moment. For ^^Cr, 
with fi = 6/iB (/^B the Bohr magneton), one has = 91 ao, E^ = 1-7 x 10~^ eV, and 
uJd/271 = 4.2 X 10^ Hz. 

In these "atomic" units, the Hartree equation of the ground state of a system of 
identical bosons in an external trapping potential VJ;rap = m{u'^r^ + u;^z^)/2, all in the 
same single-particle orbital ip, interacting via the contact interaction and the magnetic 
dipole-dipole interaction, assumes the dimensionless form 




(1) 



A + + -flz^ + N87i—\ip{r)\^ 
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+ N |^(rOP I dh'] V^(r)=£^(r). (2) 

J |r — r -I 

Here, e is the chemical potential, a is the scattering length, and 7^,^ = ujr^z/ (2t^d) are 
the dimensionless trap frequencies. Setting \1/ = \/Ntjj, one recovers the familiar form 
of the time-independent Gross-Pitaevskii equation for dipolar quantum gases. 

In the above "atomic" units, the mean-field Hamiltonian in ([2]) obeys a scaling law 
with respect to the number N of atoms: Let ipi^) be a solution of the (formal) one-boson 
problem for a given scattering length a /ad and trap frequencies jr,z, 

H^i{N = l,a/ad,7r,.)(f) V^(f) = e , (3) 

then ipi'^) •= N'^^"^ tp{i), with r = A^f, solves the A^-boson problem for the same 
scattering length a/ ad, 

H^f{N, a/ ad, jr,z) (r) ^ [r) = e (r) , (4) 

but with trap frequencies ■jr,z = Jr,z/N'^ and chemical potential e = e/N"^. Note that the 
particle number scaling leaves the aspect ratio A = jz/lr invariant. Thus, the physical 
properties of Bose condensates of dipolar quantum gases quite generally only depend on 
the value of the scattering length a/ ad and the particle number scaled trap frequencies 
N'^lr,z or, alternatively, on the aspect ratio A and the scaled geometric mean of the trap 
frequencies iV^7 = N'^'^'^^^'j^^ . We note that the dimensionless parameter D introduced 
by Dutta and Meystre [27] and Ronen et at [28] to measure the effective strength of 
the dipole interaction in trap frequency units is related to our scaling parameters by 
D = {N^-fr/^y/^. 

As an application of the particle number scaling law we emphasise that the 
experimental results reported by Koch et al. [20] for the stabilisation of dipolar chromium 
quantum gases for particle numbers ~ 20.000 and a trap frequency u/27i = 720 Hz 
correspond to a value of the scaled trap frequency A^^7 = 3.4 x 10^. Therefore they 
directly carry over to any pairs of particle numbers and mean trap frequencies with this 
value of the scaled trap frequency and the same aspect ratios! 

To study the nonlinearity effects of the time-independent Gross-Pitaevskii equation 
for dipolar gases we adopt the familiar variational ansatz of a (normalised) Gaussian 
type orbital (e. g. [9,18,20]) 

ip{r) = Aexp [-{Arr"^ + A-,z^")\ (5) 

and exploit the time- dependent variational principle for the mean-field energy to 
determine the width parameters A^ and A^. It is well known that solutions can be found 
only in certain parameter ranges and that at critical values the condensate collapses. 
What - to the best of our knowledge - for dipolar quantum gases has gone unnoticed 
before is that at these stability thresholds actually two solutions of the extended Gross- 
Pitaevskii equation ([2]) disappear. The situation is depicted in figure [1] where the 
chemical potential is plotted as a function of the scattering length for different values 
of the trap aspect ratio A and the fixed value of the scaled geometric mean of the trap 
frequency of A^^7 = 3.4 x 10^. It can be seen that, as the scattering length is increased. 
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Figure 1. Chemical potential for the mean trap frequency N'^j = 3.4 x 10'', used 
in the experiments of Koch et al. [20], and different values of the aspect ratio. The 
tangential character of the bifurcations is particularly evident for the trap aspect ratios 
A < 1.2. 




Figure 2. Stability thresholds, i.e. the critical scattering lengths Ocrit/ad, as a function 
of the scaled trap parameter 7V^7 and the aspect ratio A. In the limit N'^'j — > oo, A — > 0, 
one has acrit/od = 1/6. 

for every value of A two solutions are born in a tangent bifurcation, one corresponding 
to the ground state and the other to a collectively excited state. An inspection of 
the mean-field energies shows that the excited state corresponds to the branch of the 
chemical potential which diverges for a/a^ — > 1/6. 

The behaviour of the stability thresholds over a wide region of the parameter space 
spanned by A^^7 and A is shown in figure [2j Large values of N'^'y correspond to the regime 
where the dipole-dipole interaction is dominant, as in the experiments of Koch et al. [20]. 
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For A^^7 < 1 the contact interaction prevails, and the dipole-dipole interaction is only 
a small perturbation. 

In what follows we shall investigate the structure of the bifurcation, the stability and 
the dynamics of dipolar condensates. The first aspect we want to highlight is that dipolar 
Bose-Einstein condensates near the stability thresholds are experimental realizations of 
physical systems with "exceptional points", a phenomenon hitherto observed only in 
open quantum systems, described by non-Hermitian Hamiltonians ( [29,30], see also [31] 
and references therein). Exceptional points are positions in the parameter space where 
both the energy eigenvalues and the wave functions of (usually) two eigenstates pass 
through a branch point singularity as functions of the parameters and, consequently, 
are identical at the critical set of parameters. While in linear Schrodinger equations 
two eigenstates can only become identical if the Hamiltonian is non-Hermitian, two 
coalescing eigenstates can also occur for a nonlinear Schrodinger equation as an effect 
of its nonlinear ity. 

This indeed the case for the stationary solutions of the Gross-Pitaevskii equation 
with dipole-dipole interaction. At the stability thresholds the energy eigenvalues and the 
corresponding wave functions of the two bifurcating states pass through a branch point 
singularity of the Hamiltonian, the energies N'^e and the corresponding wave functions 
are identical. The way to reveal the branch point singularity structure is to continue 
the scattering length a/a^ into the complex plane and to check a well-known property 
of exceptional points: If in traversing a full circle around the critical value acrit/aa at 
the point of bifurcation, a/ ad = (ocrit/o-d) + fpe**^, ip = . . . 27i, a permutation of the two 
solutions occurs, an exceptional point is located within the circular area [29]. 

For the parameters of the Koch et al. [20] experiment we have performed such 
an analysis, and for the case of an almost purely dipolar quantum gas the results are 
shown in figure [3l As one goes around a small circle in the complex plane with the 
critical scattering length at the centre, the two solutions are permuted, i.e. we have 
an exceptional point. It must be said, however, that the usual way of experimentally 
proving the occurrence of an exceptional point in open quantum systems, changing two 
real physical parameters to traverse a circle in the complex energy plane, cannot be 
applied here. In the case of the Gross-Pitaevskii equation it is by continuing one real 
parameter, the scattering length, into the complex plane that the nature of the stability 
thresholds as exceptional points is revealed. 

We now analyse the dynamics of the condensate wave functions. To do so we start 
from the time- dependent Gross-Pitaevskii equation, which is ([2]) with the replacement 
e — > i^, and generalise the ansatz ([5]) to a time- dependent Gaussian type orbital 
with complex width parameters Ar{t) and A^(t) and complex normalisation factor A{t) 
(cf. [21]). We can then apply the time-dependent variational principle || iip{t)—Hip{t) |p 
= min. [32,33] to derive a system of ordinary nonlinear differential equations for the time 
evolution of the real and imaginary parts of the variational parameters Ar{t) and Az{t). 
These can be shown to be equivalent to canonical equations of motion belonging to a 
two-dimensional nonintegrable autonomous Hamiltonian system. For the time evolution 
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Figure 3. Chemical potentials N'^e (b) of the two stationary solutions emerging in the 
tangent bifurcation for a circle (a) with radius g ~ 10^^ around Ocrit/ad — —0.01929 in 
the complex extended parameter plane (iV^7 = 3.4 x lO"', A = 6.0). The permutation 
of the chemical potentials after a full circle proves the existence of an exceptional point. 



one can therefore expect all the features familiar from nonlinear dynamics studies of such 
systems, including a transition to chaos. For brevity, the details of these calculations 
are relegated to a subsequent paper. Here we shall concentrate on the results. 

We first investigate the stability of the two independent solutions of the stationary 
Gross-Pitaevskii equation which emerge from the bifurcations. To this end the four 
equations of motion for the real and imaginary parts of and are linearised around 
the values corresponding to the stationary states. In this way for each of the two states 
we obtain four eigensolutions ~ e^^ of the linearised system of equations, with 

eigenvalues k. For the ground state all eigenvalues turn out to be purely imaginary, 
proving that the state is indeed dynamically stable. For the collectively excited state 
one also finds a positive real eigenvalue, and hence the state is dynamically unstable. 
Similar behaviour was found by Huepe et al. [3, 4] in their study of the stability of 
bifurcating solutions with only an attractive contact interaction present. 

The system of nonlinear first-order differential equations also serves to investigate 
the time evolution of any initial state of the condensate by following the corresponding 
trajectories in the four- dimensional configuration space spanned by the coordinates of 
the real and imaginary parts of A^ and A-^. Since the total mean-field energy is a constant 
of motion the trajectories are restricted to three-dimensional hyperplanes, and their 
behaviour can most conveniently be visualised by two-dimensional Poincare surfaces of 
section defined by requiring one of the coordinates to assume a fixed value. 

We consider Poincare surfaces of section defined by the condition that the imaginary 
part of is zero. Each time the trajectory crosses the plane Im(y42) = 0, the real and 
imaginary parts of Ar(^) = ^r(^) +^^r(^) recorded. In figure H] surfaces of section are 
plotted for seven different, increasing, values of the mean-field energy. Again the physical 
parameters of the experiment of Koch et al. [20] are adopted, and the scattering length 
is fixed to aja^ = 0.1, away from its critical value. At these parameters, the variational 
mean-field energy of the ground state is NEgs = 4.24 x 10^ and represents the local 
minimum on the two-dimensional mean-field energy landscape, plotted as a function of 
the (real) width parameters. The variational energy of the second, unstable, stationary 
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Figure 4. Poincare surfaces of section of the condensate wave functions represented 
by tlieir widtli parameters for seven different mean-field energies at tlie scaled trap 
frequency N'^^ — 3.4 x 10'*, aspect ratio A = 6, and the scattering length a/ad — 0.1. 
The surfaces of section labelled 1 to 7 correspond to the increasing values of the mean- 
field energy NE = 4.24 x 10^5.00 x 10^ 6.00 x 10^9.00 x 10^1.50 x 10^,3.00 x 
10*^, 6.00 X 10*^, respectively. 

state at these experimental parameters is NE^s = 6.24 x 10^, it corresponds to the saddle 
point on the mean-field energy surface. Between these two energy values the motion on 
the trajectories is bound, while for energies above the saddle point energy the motion on 
the trajectories can become unbound: once the saddle point is traversed by a trajectory 
Ar{t), Az{t), the parameters run to infinity, Al.{t),Al{t) —>■ oo, meaning a shrinking of 
the quantum state to vanishing width, i.e. a collapse of the condensate takes place. 

The surface of section labelled 1 in figure H] belongs to the energy of the ground state. 
At this energy the kinematically allowed region for the crossing points of the trajectories 
is confined to a single stable stationary point. At the next higher energy (surface of 
section labelled 2) the kinematically accessible region in configuration space has grown, 
and the initially stationary state has evolved into a periodic orbit (fixed point in the 
surface of section), corresponding to a state of the condensate whose motion is periodic. 
The oscillations of the width parameters Ar{t) and A^it) represent oscillatory stretchings 
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of the condensate along the r and z directions. The stable periodic orbit in the surface 
of section is surrounded by elliptical, quasi-periodic orbits, representing quasi-periodic 
oscillations of the condensate. The surface of section 3, at the next higher energy, 
reveals that bifurcations have occurred, creating new stable and unstable periodic states, 
manifested by the emergence of additional elliptical islands and separatrices in the 
surface of section. 

The surface of section labelled 4 is the first in figure H] with an energy value above 
the saddle point energy. Now chaotic orbits have appeared which surround the stable 
regions. In contrast to the (quasi-) periodic stretching oscillations of the condensate 
within the elliptical islands, the chaotic motion of the parameters describes a condensate 
which does not yet collapse but whose widths fluctuate irregularly. 

One might imagine that well above the saddle point energy stable condensate wave 
functions no longer exist. However, in the surfaces of section labelled 5, 6, and 7 
regular islands are still clearly visible. These stable islands are surrounded by chaotic 
trajectories. Since ergodic motion along these trajectories comes close to every point 
in the configuration space, the chaotic motion sooner or later leads to a crossing of the 
saddle point and then to the collapse of the condensate wave functions. It can be seen 
that with growing energy above the saddle point the sizes of the stable regions gradually 
shrink. The reason why the kinematically allowed regions surrounding the stable islands 
are hardly recognisable any more in these surfaces of section is that high above the saddle 
point energy the chaotic motion becomes more and more unbound, which means that 
the trajectories cross the Poincare surfaces of section only a few times, if ever, before 
they escape to infinity and collapse takes place. 

It must be emphasised, however, that stable islands do persist even far above the 
saddle point energy, implying the existence of quasi-periodically oscillating nondecaying 
modes of the condensate wave functions. 

The prediction of such modes of energetically excited solutions of the time- 
dependent Gross-Pitaevskii equation for cold dipolar quantum gases is a result of our 
analysis. It would certainly be an intriguing and challenging task to examine whether 
it is possible to prepare excited states of dipolar quantum gases of this type in both 
the regular as well as in the chaotic regions, to distinguish between the two different 
dynamics, and to access the stable regions high above the energy of the unstable 
stationary state. One way of creating the collectively excited states one might imagine 
is to prepare the condensate in the ground state, and then to non-adiabatically reduce 
the trap frequencies. Clearly experimental investigations along these lines are strongly 
encouraged. 

One might ask, however, whether the Gross-Pitaevskii equation underlying our 
calculations is adequate at all to describe complex dynamics of this type in real dipolar 
quantum gases. In particular in the chaotic regime local density maxima might occur for 
which losses by two-body or three-body collisions would have to be taken into account. 
However, by virtue of the scaling law (jl]) parameter ranges can always be found where the 
particle densities remain small even in these regimes and the Gross-Pitaevskii equation 
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is applicable. 

In this paper we have used a simple variational ansatz ([5]) to model the dynamics of 
axisymmetric solutions of the Gross-Pitaevskii equation of dipolar gases. The advantage 
of the variational technique is that the analysis of the nonlinearity effects becomes 
especially transparent. Of course the results must be checked by accurate numerical 
simulations. Such simulations, and the extensions to fully three-dimensional and 
structured condensate wave functions [27,34], are under way. Comparisons of variational 
and accurate numerical results [24-26] in an alternative system with a long-range 
interaction, viz. Bose condensates with an attractive gravity-like 1/r interaction [23], 
have shown that the nonlinear dynamical properties found in the variational calculation 
are recovered in the numerical calculations, and that the quantum behaviour may be 
even richer. As a common feature of the propagation of quantum wave functions, the 
real-time evolution of arbitrary condensate wave functions of dipolar gases will exhibit 
complicated fluctuations. An interpretation of their full dynamics, and in particular the 
search for periodic, quasiperiodic or chaotic structures, therefore will hardly be possible 
without the guidance of the variational results obtained in this paper. 
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